PrecipitationDailyRead Subroutine

SUBROUTINE PrecipitationDailyRead &
!
( time, dem)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time !!current time
TYPE (grid_real), INTENT(IN) :: dem

!local declarations:
TYPE (DateTime) :: time_toread !! time to start reading from
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname
INTEGER (KIND = short) :: i, j

!-------------------------end of declarations----------------------------------

IF ( .NOT. (time < timeNewPrecipitation) ) THEN
   
   !time_toread = time + - (dtPrecipitationDaily - raingagesDaily % timeIncrement)
   time_toread = time
   timeString = time
   CALL Catch ('info', 'PrecipitationDaily',   &
		                 'read new precipitation data: ', &
                     argument = timeString)
   
   !update lapse rate when required
    IF (elevationDrift == 1) THEN
      IF (time == lapse_map % next_time) THEN !update lapse_map
          varname = lapse_map %var_name
          filename = lapse_map % file_name
          !destroy current grid
          CALL GridDestroy (lapse_map)
          !read new grid in netcdf file
          CALL NewGrid (layer = lapse_map, fileName = TRIM(filename), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(varname), &
                         time = time_toread)
      END IF
    END IF

   SELECT CASE (interpolationMethod)
    CASE (0) !input grid
      
		CALL ReadField (fileGrid,  time_toread, &
        dtPrecipitationDaily, dtGrid, &
        'C', precipitationRateDaily, &
            varName = precipitationRateDaily % var_name)
         
    CASE DEFAULT !requires interpolation
      !read new precipitation data
      CALL ReadData (network = raingagesDaily, fileunit = fileunit, &
                      time = time_toread, aggr_time = dtPrecipitationDaily, &
                      aggr_type = 'sum', tresh = valid_prcn)
      
      IF (elevationDrift == 1) THEN

        !shift precipitation observation to reference elevation by applying lapse rate
        CALL ShiftMeteoWithLapse (raingagesDaily, lapse_map, refElevation, &
                                  stationsRefElev, dtPrecipitationDaily)
        
        !interpolate 
        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = precipitationRateDaily, &
                            devst = grid_devst)                  
                
        ELSE 
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = interpolatedMap (i), &
                            devst = grid_devst) 
                      
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        precipitationRateDaily % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO
        END IF

        !Shift back interpolated field to terrain surface
        CALL ShiftBackWithLapse (precipitationRateDaily, dem, &
                                lapse_map, refElevation, &
                                dtPrecipitationDaily)
         
      ELSE
            IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                 CALL Interpolate (network = raingagesDaily, &
                              method = interpolationMethod, &
                              near = neighbors, &
                              idw_power = idw_power, & 
                              anisotropy = krige_anisotropy, &
                              varmodel = krige_varmodel, &
                              lags = krige_lags, &
                              maxlag = krige_maxlag, &
                              grid = precipitationRateDaily, &
                              devst = grid_devst) 
                
            ELSE
                !loop trough interpolation methods
               DO i = 1, 3
                  IF (interpolationMethod_vector(i) > 0) THEN
                      
                       CALL Interpolate (network = raingagesDaily, &
                              method = interpolationMethod_vector(i), &
                              near = neighbors, &
                              idw_power = idw_power, & 
                              anisotropy = krige_anisotropy, &
                              varmodel = krige_varmodel, &
                              lags = krige_lags, &
                              maxlag = krige_maxlag, &
                              grid =  interpolatedMap (i), &
                              devst = grid_devst) 
                    
                  END IF
               END DO

               !compose the final interpolated field
               DO i = 1, interpolationMethod_map % idim
                  DO j = 1, interpolationMethod_map % jdim
                      IF (interpolationMethod_map % mat(i,j) /= &
                          interpolationMethod_map % nodata ) THEN
                         precipitationRateDaily % mat (i,j) = &
                         interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                      END IF
                  END DO
               END DO

            END IF  !1 or more interpolation methods       
      END IF  !elevationDrift
		
	END SELECT

  !apply scale factor and offset, if defined
	IF (offset_value /= MISSING_DEF_REAL) THEN
	   CALL Offset (precipitationRateDaily, offset_value)
	END IF
	
	
	IF (scale_factor /= MISSING_DEF_REAL) THEN
	   CALL Scale (precipitationRateDaily, scale_factor)
    END IF

   
  !filter data < 0.
    DO i = 1, precipitationRateDaily % idim
      DO j = 1, precipitationRateDaily % jdim
           IF (precipitationRateDaily % mat(i,j) /= precipitationRateDaily % nodata) THEN
              IF ( precipitationRateDaily % mat(i,j) < 0.) precipitationRateDaily % mat(i,j) = 0.
           END IF
      END DO
    END DO
    
  !grid exporting
  IF(export > 0 ) THEN
	  IF( time == timeNewExport .AND. &
        time >= export_start .AND. &
        time <= export_stop ) THEN
        timeString = time
        timeString = timeString(1:19)
        timeString(14:14) = '-'
		    timeString(17:17) = '-'
        
        grid_devst % reference_time = precipitationRateDaily % reference_time
        IF (needConversion) THEN
           CALL GridConvert (precipitationRateDaily, exportedGrid)
           CALL GridConvert (grid_devst, exportedGridVar)
        ELSE
           exportedGrid = precipitationRateDaily
           exportedGridVar = grid_devst
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_precipitation.asc'
              CALL Catch ('info', 'PrecipitationDaily',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII)
              
              IF (krige_var == 1) THEN !export kriging variance
                    export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                           '_precipitation_variance.asc'
                    CALL Catch ('info', 'PrecipitationDaily',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                    CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII)
              END IF
              
        CASE (2) !esri-binary
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_precipitation'
              CALL Catch ('info', 'PrecipitationDaily',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY)
              
              IF (krige_var == 1) THEN !export kriging variance
                   export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                               '_precipitation_variance'
                   CALL Catch ('info', 'PrecipitationDaily',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                   CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY)
              END IF
              
        CASE (3) !net_cdf
              CALL SetCurrentTime (time, exportedGrid)
              CALL Catch ('info', 'PrecipitationDaily',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, 'precipitation_amount', 'append')
              
              IF (krige_var == 1) THEN !export kriging variance
                  CALL SetCurrentTime (time, exportedGridVar)
                  CALL Catch ('info', 'PrecipitationDaily',   &
		                       'exporting grid to file: ', &
                              argument = TRIM(export_file_var))
                 CALL ExportGrid (exportedGridVar, export_file_var, 'precipitation_amount_variance', 'append')
              END IF
        END SELECT
        timeNewExport = timeNewExport + export_dt
    END IF
  ENDIF


   
  !conversion mm => m/s
  DO i = 1, precipitationRateDaily % idim
    DO j = 1, precipitationRateDaily % jdim
       IF (precipitationRateDaily % mat(i,j) /= precipitationRateDaily % nodata) THEN
           precipitationRateDaily % mat(i,j) = precipitationRateDaily % mat(i,j) / &
           dtPrecipitationDaily / 1000.
       END IF
    END DO
  END DO


  !time forward
  timeNewPrecipitation = timeNewPrecipitation + dtPrecipitationDaily
END IF

RETURN
END SUBROUTINE PrecipitationDailyRead